* Program to construct bootstrap percentiles of noise-signal ratio from real-time gap
* data over the 1966-1997 sample period.
* Note:  blockwidth and nobs are set by hand, and this program requires the number
* of blocks to divide evenly into the sample for the circular MBB procedure.

* THIS VERSION:  Uses block width of 5 (and 26 blocks), so bootstrap sample has
*                130 observations (note that this requires hardwiring certain parts
*                of the code that use NOBS to equal 130 instead--that is, NOBS is
*                now the size of the bootstrap sample, and 128 is the number of
*                entries in the original dataset).

* This program:  Uses a block width of 5 and the RMSE-based noise-signal ratio


calendar 1947 1 4
allocate 2013:04+5000

seed 4892

open data table1_boot_data.rat
data(format=rats) / gb_gap_rt gb_gap_gdp99q3

* GB_GAP_RT:       Real-time gap estimate
* GB_GAP_GDP99Q3:  Final gap for 1966-1997 period


* Put the relevant data series into a vector (note that
* the make command will actually yield an array...)

make rtgap_vec 1966:01 1997:04
# gb_gap_rt
disp %rows(rtgap_vec); * This will be 128, and will no longer equal NOBS, below

make fingap_vec 1966:01 1997:04
# gb_gap_gdp99q3
disp %rows(fingap_vec); * This will be 128, and will no longer equal NOBS, below

compute gaprev = fingap_vec-rtgap_vec
compute sqgaprev = gaprev.^2.0

compute blockwidth=5
compute nobs=130
compute numblks=26
compute numdraws=5000

declare vector xxx(numblks)
declare vector[integer] entries(nobs)
declare vector b_sqrevs(nobs)
declare vector b_fin(nobs)
declare vector n_s_ratios(numdraws)



do drawiter=1,numdraws


ewise xxx(i)=%raninteger(1,128)

* Use these two lines of code (with NUMDRAWS=1) to test the calculations (they just return the original series)
* compute qqq=%seqa(1,blockwidth,nobs)
* ewise xxx(i)=qqq(i)


do iter=1,numblks
  do subiter=1,blockwidth
    compute entries(((iter-1)*blockwidth)+subiter)=%clock(fix(xxx(iter))+(subiter-1),128)
  end do subiter
end do iter

do iter=1,nobs
  compute b_sqrevs(iter)=sqgaprev(entries(iter),1)
  compute b_fin(iter)=fingap_vec(entries(iter),1)
end do iter

compute rmse=sqrt(%avg(b_sqrevs))
compute sqdev=(b_fin-%avg(b_fin)).^2.0
compute nnn=%rows(sqdev)
compute nscale=float(nnn)/float(nnn-1)
compute stddev=sqrt(nscale*%avg(sqdev))

compute n_s_ratios(drawiter)=rmse/stddev

end do drawiter



set nsr_series 1 numdraws = n_s_ratios(t)
statistics(fract) nsr_series

